############################################################
############################################################
##### FIGURE REPLICATION CODE FOR: 
##### IN THE SHADOWS OF GREAT MEN
##### RETIRED LEADERS AND INFORMAL POWER CONSTRAINTS IN AUTOCRACIES
############################################################
############################################################


library(ggplot2)
library(Cairo)
library(foreign)
library(gridExtra)
library(gtable)
library(scales)
library(data.table)
library(tidyverse)
library(openxlsx)
library(haven)
library(readstata13)
Sys.setlocale(,"CHS")


setwd("D:/Dropbox/Academic/Research/World leader power/Replication")
listdele<-function(x){
  a<-x[complete.cases(x),]
  return(a)
}

theme_set(theme_bw() + theme(legend.position="bottom",legend.direction="horizontal",
                             axis.title.x=element_text(vjust=-.5),
                             axis.text.y =  element_text(colour = "black", hjust = 0 , vjust=.5 ),
                             title=element_text(face="bold",vjust=1.5),
                             strip.text=element_text(size=13),
                             plot.title = element_text(hjust = 0.5)))


pd<-position_dodge(width=0.4)
cl<-c("red","steelblue")

source("stat_binscatter.r")

d<-read_dta("replication_data_country_year.dta")


################################################################################
#### Figure 1. Comparing Ngram-based Power Index with Existing Measures ########
################################################################################

# panel 1 #
gh<-ggplot(d[is.na(d$t1_personalism)==F,],aes(x=rpbook_noc))+geom_histogram()+xlab("")+ylab("")+ggtitle("1. Power Index (Ngram)")

# panel 2 #
ggwf<-ggplot(d,aes(x=t1_personalism))+geom_histogram()+xlab("")+ylab("")+ggtitle("2. Personalism Index (GWF)")

# panel 3 #
gggs<-ggplot(d,aes(x=xhatmean))+geom_histogram()+xlab("")+ylab("")+ggtitle("3. Power Consolidation Index (GS)")


# panel 4. Ngram vs. GWF #
g.comp2<-ggplot(d, aes(x=rpbook_noc   , y = t1_personalism)) +
  stat_binscatter(bins = 10,geom="pointrange")+
  xlab("Power index (Ngram)")+ylab("Personalism (GWF)")+
  ggtitle("4. Ngram vs. GWF")+
  annotate("text",x=-2.1,y=0.59,label='beta == "0.273***"',parse=T)

# panel 5. Ngram vs. GS #
g.comp3<-ggplot(d, aes(x= rpbook_noc   , y= xhatmean)) +
  stat_binscatter(bins = 10,geom="pointrange")+
  xlab("Power index (Ngram)")+ylab("Power consolidation (GS)")+
  ggtitle("5. Ngram vs. GS")+
  annotate("text",x=-2.1,y=0.59,label="beta == 0.043",parse=T)


# panel 6. by polity score #
d$type<-NA
d$type[which(d$polity2<=-6)]<-1
d$type[which(d$polity2>-6 & d$polity2<=0)]<-2
d$type[which(d$polity2>0 & d$polity2<=5)]<-3
d$type[which(d$polity2>6)]<-4

dp<-rbindlist(lapply(1:4,function(x) mean_se(d$rpbook_noc[which(d$type==x & d$year<=2010)],mult = 1.96)))
dp$type<-factor(1:4,levels=1:4,labels=c("Autocracy\n(-10 to -6)",
                                        "Hybrid regime \n(-5 to 0)",
                                        "Hybrid regime \n(1 to 5)",
                                        "Democracy\n(6 to 10)"))

gp<-ggplot(dp,aes(x=type,y=y))+geom_pointrange(aes(ymin=ymin,ymax=ymax),size=.4)+ylab("Power index (Ngram)")+xlab("Polity")+
  ggtitle("6. Ngram by Polity Score")+
  scale_y_continuous(limits=c(-1.2,1.2))

# panel 7. by regime type #

personal<-mean_se(d$rpbook_noc[which(d$gwf_personal2==1 & d$auto_samp==1)],mult = 1.96)

party<-mean_se(d$rpbook_noc[which(d$gwf_party2==1& d$auto_samp==1)],mult = 1.96)

military<-mean_se(d$rpbook_noc[which(d$gwf_military2==1& d$auto_samp==1)],mult = 1.96)

monarchy<-mean_se(d$rpbook_noc[which(d$gwf_monarchy2==1& d$auto_samp==1)],mult = 1.96)


dreg<-rbind(personal,party,military,monarchy)
dreg$regime<-factor(c("Personalist\n","Party\n","Military\n","Monarchy\n"))
dreg<-dreg[order(dreg$y),]
dreg$regime<-factor(dreg$regime,levels=dreg$regime)


gregim<-ggplot(dreg,aes(x=regime,y=y))+geom_pointrange(aes(ymin=ymin,ymax=ymax),size=.4)+ylab("Power index (Ngram)")+xlab("GWF Regime Type")+  
  ggtitle("7. Ngram by Autocratic Regime Type")+
  scale_y_continuous(limits=c(-1.2,1.3))




# combine everythign #
pdf("intext_validation.pdf",height=9,width=9)
grid.arrange(gh,ggwf,gggs,g.comp2,g.comp3,gp,gregim,ncol=2,nrow=4,layout_matrix =rbind(c(1,1),c(2,3),c(4,5),c(6,7)))
dev.off()


png("intext_validation.png",height=8000,width=8000,res = 1000)
grid.arrange(gh,ggwf,gggs,g.comp2,g.comp3,gp,gregim,ncol=2,nrow=4,layout_matrix =rbind(c(1,1),c(2,3),c(4,5),c(6,7)))
dev.off()



####################################################################
#############  Figure 2. Variation in Predecessor Power  ###########
####################################################################

d$countryn[which(d$countryn=="malaysia")]<-"Malaysia"
d$countryn[which(d$countryn=="eygpt")]<-"Egypt"
d$cenX<-d$cen
d$cenX[which(d$countryn=="Russia" & d$year>=2000 & d$year<=2007)]<-"Putin-A"
d$cenX[which(d$cenX=="Lee Hsien Loong")]<-"LEE Hsien Loong"

ct<-c("China","Mexico","Vietnam","Tanzania","Russia","Singapore","Malaysia","Cuba" )


dc<-d[which(d$countryn%in%ct),]
dc<-dc[-which((dc$main_samp==0 & dc$countryn%in%c("China","Vietnam","Tanzania","Egypt","Singapore","Malaysia","Cuba"))|dc$cen=="Boris Yeltsin"),]

dc<-dc[-which((dc$year>1999 & dc$countryn=="Mexico")|(dc$countryn=="Malaysia"&dc$year>=2018)|(dc$year%in%c(2011,2012,2013) & dc$countryn=="Egypt")|(dc$countryn=="Senegal" & (dc$year<1960|dc$year>=2000))),]

ds<-dc %>% group_by(cenX) %>%
  mutate( xmin=min(year),
          xmax=max(year)) %>%
  arrange(cowcode, year) %>%
  mutate(period = factor(paste(xmin,xmax,countryn,sep="-")),
         fill = NA)


for(i in 1:length(ct)){
  ds$fill[which(ds$countryn==ct[i])]<-as.numeric(factor(paste(ds$xmin[which(ds$countryn==ct[i])],ds$xmax[which(ds$countryn==ct[i])],sep="-")))%%2
}


ds[which(ds$year!=ds$xmin),c("xmin","xmax")]<-NA

gl<-ggplot(ds,aes(x=year,group=period))+
  geom_line(aes(y=log_maxOfallprepower,color="1",linetype="1"))+
  geom_line(aes(x=year,y=rpbook_noc,color="2",linetype="2"))+
  geom_rect(aes(ymin=-Inf,ymax=Inf,xmin=xmin,xmax=xmax+1,fill=factor(fill)),alpha=.3)+facet_wrap(.~countryn,ncol=2,scale="free_y")+
  scale_fill_manual("",values=c("steelblue","orange"))+
  guides(fill="none")+xlab("")+ylab("Incumbent's and predecessor's power")+
  scale_color_manual("",values=c("black","red"),labels=c("Predecessor's power","Incumbent's power"))+
  scale_linetype_manual("",values=c("dashed","solid"),labels=c("Predecessor's power","Incumbent's power")) 



pdf("id_variation.pdf",height=5.6,width=8)
gl
dev.off()

png("id_variation.png",height=5600,width=8000,res=1000)
gl
dev.off()

###################################################################
################### Figure 3. Parallel Trends #####################
###################################################################


dt<-listdele(read.csv("parallel.csv",stringsAsFactors=F))
dt$label<-c("In 3\n years",
            "In 2\n years",
            "In 1\n year",
            "1 year ago",
            "2 years ago",
            "3 years ago")

rect<-data.frame(xstart=c(0,3.5),xend=c(3.5,7),col=letters[1:2])
rect$x.text<-ifelse(is.infinite(rect$xend),9,rect$xend)
gtext<-data.frame(txt=c("Will Lose a Predecessor","Predecessor Lost"),x=(rect$xstart+rect$x.text)/2,y=0.5)
dt$xax<-factor(dt$label,levels=dt$label)

gpara<-ggplot()+
  geom_rect(data=rect,aes(xmin=xstart,xmax=xend,ymin=-Inf,ymax=Inf,fill=col),alpha=.2)+
  geom_point(data=dt,aes(x=xax,y=beta),size=2)+
  geom_errorbar(data=dt,aes(x=xax,ymin=lb90,ymax=ub90,y=beta),width=.1)+
  geom_hline(yintercept=0,linetype="dashed",color="red")+
  scale_x_discrete(labels=dt$label)+
  scale_fill_manual(values=c("grey60","white","grey60"))+
  xlab("")+ylab("Estimated Effect")+
  geom_text(data=gtext,aes(label=txt,x=x,y=y),size=5)+
  guides(fill="none")


pdf("parallel.pdf",width=7.5,height=4)
gpara
dev.off()

png("parallel.png",width=7500,height=4000,res=1000)
gpara
dev.off()

####################################################################
############## Figure 4. Results bu Regime Type  ###################
####################################################################


#### Only showing results by regime type #############


ds<-listdele(read.csv("sub_coef.csv",stringsAsFactors = F)) 
ds$lab<-c("Party-based regimes\n(n = 779)","Military regimes\n(n = 180)","Personalist regimes\n(n = 220)")

ds$coef<-paste(round(ds$beta,2),"\n"," [",round(ds$lb90,2),",",round(ds$ub90,2),"]",sep="")


ds$lab<-factor(ds$lab,levels=ds$lab)


g2<-ggplot(ds,aes(x=lab,y=beta))+
  geom_pointrange(aes(ymin=lb90,ymax=ub90))+
  geom_hline(yintercept=0, linetype="dashed",colour="red")+
  xlab("")+ylab("Estimated effect")+
  scale_y_continuous(limits=c(-1.5,1.2))+
  geom_text(aes(label=coef,x=as.numeric(lab)+.24))+
  theme(axis.text.x.bottom = element_text(size=12))


pdf("by_regime.pdf",height=3,width=6.5)
g2
dev.off()


png("by_regime.png",height=3000,width=6500,res=1000)
g2
dev.off()

